Melody: a note detection system

This is a program which tries to recognize notes from a melody and produce a visualization of it. Here is an example of what the result looks like with a short clip from a Beethoven sonata (the audio can be found below in the "Note onset detection" section of this notebook):

There is a section at the bottom of the notebook with more examples. "Run all cells" may be helpful here since the examples require all the previous code cells to be executed first.

The general strategy

We're using a Fourier Transform based strategy where we break the audio into small segments and look at the FT of each segment to determine which note is being played in each. Currently, our approach doesn't generalize well to polyphonic music and is limited to pieces with one note being played at a time.

Before running on actual music, we "train" on a chromatic scale and save the FTs of each note in the scale. We'll use these FTs as reference points when trying to detect pitch. Going through this "training" instead of using a fixed frequency-to-note mapping can help to calibrate to the specific timbre and tuning of different instruments.

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation
import scipy.io.wavfile, scipy.fft, scipy.signal, scipy.stats
import IPython

%matplotlib inline
%config InlineBackend.figure_format = 'svg'

I. Note onset detection

We want to find the timestamps of the start of each note (this is only used for the training phase, when we break the chromatic scale down into individual notes).

As an example, there are nine notes in this clip (Beethoven Piano Sonata No. 31, mvt. 3, fugue theme), and we want to mark the onset of each of these.

sample_file = 'yamaha/LvB_fugue_theme.wav'
IPython.display.Audio(sample_file)

Our approach involves picking out peaks in the raw audio waveform. Since each note produces a new sound, this tends to correlate to a spike in the audio. We're interested in these macro-level peaks:

sample_rate, data = scipy.io.wavfile.read(sample_file)
data = np.abs(data)
plt.plot(data);

But there are many other peaks (on the "micro"-level) caused by the sinusoidal nature of sound. We can see these micro-level peaks if we zoom in on a small part of the audio:

plt.plot(data[73000:74000]);

We want to work with a "zoomed out" version the audio so that we only see the macro-level peaks and ignore the micro-level ones. To do this, we can slide a window over the waveform, take the maximum amplitude of the audio within that window, and discard the rest of the data. This achieves "zooming out" in a sense by hiding the micro-level movements in the waveform and keeps only the general, macro trend.

The condense() function performs this step, returning the condensed data (orange line in the graph below) as well as timestamps in the original audio corresponding to the maximums of each window to be used later (green dots in the graph below).

def condense(data, window_size):
    padded = np.pad(data, (0, (window_size - len(data) % window_size) % window_size))
    reshaped = padded.reshape((-1, window_size))
    return reshaped.max(axis = 1), reshaped.argmax(axis = 1) + np.arange(0, len(data), window_size)

c_window = 5500
c, i = condense(data, c_window)
plt.plot(data)
plt.plot(np.arange(0, len(data), c_window), c)
plt.plot(i, data[i], '.');

Finally, we can pick out the peaks from the orange line (keeping only ones that are above a certain threshold) and, if necessary, merge adjacent peaks that are too close together.

def get_peaks(data, threshold, condense_window):
    data, peaks = condense(np.abs(data), condense_window)
    for i in range(1, len(data) - 1):
        if (data[i-1] < data[i] > data[i+1] and data[i] > threshold):
            yield peaks[i]

def merge_peaks(peaks, merge_window):
    prev = -merge_window
    for p in peaks:
        if p >= prev + merge_window:
            yield p
        prev = p

Final result for onset detection:

thresh = 5000
plt.plot(data)
markers = list(merge_peaks(get_peaks(data, threshold = thresh, condense_window = c_window), 
                           merge_window = c_window * 2))
plt.plot(markers, data[markers], '.');

II. Finding reference FTs

We take the chromatic scale and find the Fourier Transform of each note in it, storing the results in the list notes.

The Fourier Transforms are cut off after a certain frequency above which there shouldn't be much activity for musical notes in a reasonable range.

def fourier_transform(data, cutoff):
    return scipy.fft.fft(data)[:cutoff]

def get_fts(data, pos, window_size, cutoff):
    return [fourier_transform(data[x:x+window_size], cutoff) for x in pos]

ft_window_hz = 8 # i.e. the window passed to Fourier spans 1/8 of a second
freq_max = 5000 # maximum frequency to include in FTs
training_data = 'yamaha/med.wav'

sample_rate, data = scipy.io.wavfile.read(training_data)
ft_window_size = sample_rate // ft_window_hz
freq_cutoff_index = freq_max // ft_window_hz

onsets = list(merge_peaks(get_peaks(data, threshold = thresh, condense_window = c_window), 
                          merge_window = 2 * c_window))
notes = get_fts(data, onsets, ft_window_size, freq_cutoff_index)
freqs = scipy.fft.fftfreq(ft_window_size, d = 1 / sample_rate)[:freq_cutoff_index]

The resulting FT for the first note in the scale (C3):

plt.plot(freqs, np.abs(notes[0]));

We're mostly interested in the peaks in the above graph (representing the fundamental frequency of 131 Hz and its overtones), so we want to produce a "sharpened" version of this FT containing only these peaks.

We can reuse the peak-picking function get_peaks() used in onset detection:

def normalize(data):
    return data / np.max(data)

def ft_peaks(data, thresh = 0.02, window = 5):
    return list(merge_peaks(get_peaks(data, threshold = thresh, condense_window = window), 
                            merge_window = window))

def sharpen(data):
    peaks = ft_peaks(data)
    sharp = np.zeros(data.shape)
    sharp[peaks] = data[peaks]
    return sharp

ref_fts = [sharpen(normalize(np.abs(ft))) for ft in notes]
ref_mags = np.sum(np.square(ref_fts), axis = 1)

Result after sharpening (of the same note C3):

plt.plot(freqs, ref_fts[0]);

III. Putting everything to use: note detection

For the given audio clip, we can slide a small window over it and find the FT of each window. We then compare this FT with the ref_fts to find the most likely note. We're treating the FT of this window as a vector $\mathbf{v}$ and finding the reference FT $\mathbf{r}$ that matches its direction the most (maximizing $\frac{\mathbf{v} \cdot \mathbf{r}}{\| \mathbf{r} \|}$).

The Fourier Transform is a linear transformation that expresses a signal in the basis of sinusoids of different frequencies; this note detection step is also based on a linear change-of-basis which switches from a basis of frequency to a basis of musical notes. Previously, we have tried to combine these two transformations, skipping the Fourier step and directly projecting the signal onto vectors representing raw waveforms of each note. However, the results were not nearly as promising. Part of this difference may likely be due to the FT sharpening step we perform between the two transformations.

def get_note(references, magnitudes, ft, threshold = 3000000):
    projected = np.dot(references, ft) / magnitudes
    note = np.argmax(projected)
    if projected[note] < threshold: return -1;
    return octave_correction(references, ft, note)

Notes that are an octave apart are harder to distinguish, so we add an "octave correction" step at the end of our detection function. This step finds the most prominent frequency in the FT (including its overtones) and compares it to the frequency in the reference FT to try to determine the right octave (each octave up doubles the frequency).

The find_strongest_freq() function estimates the frequency of the note in a given FT. For each peak in the FT, the function estimates its "strength" by adding up the strengths of all integer multiples of that frequency (with more weight given to the first few harmonics). It's roughly accurate but not quite enough to be used directly for note detection, so we only use it for octave correction.

def octave_correction(references, ft, note):
    this_freq = find_strongest_freq(ft)
    ratio = round(np.log2(this_freq / find_strongest_freq(references[note%12])))
    return note % 12 + ratio * 12

def find_strongest_freq(ft, decay = 0.2):
    maxstrength = 0
    freq = -1
    peaks = ft_peaks(ft)
    for i in range(len(peaks)):
        strength = 0
        weight = 1
        for j in range(i, len(peaks)):
            if is_almost_divisible(peaks[j], peaks[i]):
                strength += ft[peaks[j]] * weight
                weight = max(weight - decay, decay)
        if strength > maxstrength:
            maxstrength = strength
            freq = peaks[i]
    return freq

def is_almost_divisible(x, y, err_tolerance = 0.05):
    approx = y * round(x / y)
    return abs(x - approx) / approx <= err_tolerance

note_names = {'b': 'C D♭ D E♭ E F G♭ G A♭ A B♭ B'.split(), 
              '#': 'C C♯ D D♯ E F F♯ G G♯ A A♯ B'.split()}
note_colors = '''#ffb8b8 #ffc8b8 #ffdcb8 #fff1b8 #fffdb8 #dcffb8 #b8ffc0
                 #b8f8ff #b8e1ff #b8cdff #b8baff #e7b8ff #ffb8e0'''.split()

step_hz = ft_window_hz * 4
step_size = sample_rate // step_hz

The final function for running note detection on a file and producing the visualization:

def plot_notes(fname, sharp_or_flat, smallest_label = 0.3, linewidth = 2):
    sample_rate, data = scipy.io.wavfile.read(fname)
    windows = list(range(0, len(data) - step_size, step_size))
    fts = get_fts(data, windows, sample_rate // ft_window_hz, freq_cutoff_index)
    score = np.fromiter((get_note(ref_fts, ref_mags, np.abs(v)) for v in fts), int)
    segs = score2segments(score, windows)
    
    highest = max(segs, key = lambda seg: seg[2])[2]
    lowest = min(segs, key = lambda seg: seg[2])[2]

    label_threshold = len(data) // sample_rate * smallest_label
    num2name = note_names[sharp_or_flat]

    fig = plt.figure(figsize = (len(segs) // 2, (highest - lowest) // 6))
    plt.axis('off')
    for seg in segs:
        start, end, note = seg
        plt.plot([start, end], [note, note], color = note_colors[note % 12], lw = linewidth)
        if end - start > label_threshold:
            plt.annotate(num2name[note % 12], ((start+end)/2, note), va = 'bottom', ha = 'center')

def score2segments(score, timestamps, cleanup_rad = 9):
    plot_segments = []
    prev = -1
    for i in range(len(score)):
        if score[i] != -1:
            neighborhood = [score[j] for j in range(i-cleanup_rad, i+cleanup_rad+1) 
                            if j < len(score) and score[j] % 12 == score[i] % 12]
            mresult = scipy.stats.mode(neighborhood)
            note = int(mresult[0])
            if note != prev:
                plot_segments.append([timestamps[i], timestamps[i], note])
            else:
                plot_segments[-1][1] = timestamps[i]
            prev = note
        else: prev = -1;
    return plot_segments

IV. Examples

Alan Menken, "The Bells of Notre Dame" from The Hunchback of Notre Dame soundtrack:

f = 'yamaha/notre_dame.wav'
plot_notes(f, 'b')
IPython.display.Audio(f)

Jonathan Coulton, "Still Alive" from the Portal soundtrack:

f = 'yamaha/portal.wav'
plot_notes(f, '#')
IPython.display.Audio(f)

A chromatic scale (one that's faster than the one used for "training") to test a range of notes:

f = 'yamaha/fast.wav'
plot_notes(f, 'b')
IPython.display.Audio(f)

Beethoven, Piano Sonata No. 32, mvt. 2 (here, the first six notes are new notes for which we don't have reference FTs for; they are in a higher octave which was not included in the training scale, and the octave correction step was able to help correctly identify them):

f = 'yamaha/LvB_s32.wav'
plot_notes(f, '#')
IPython.display.Audio(f)

Beethoven, Symphony No. 9, mvt. 4, "Ode to Joy" theme:

f = 'yamaha/LvB_s9.wav'
plot_notes(f, '#')
IPython.display.Audio(f)

William Bolcom, "Graceful Ghost Rag":

f = 'yamaha/graceful_ghost.wav'
plot_notes(f, 'b')
IPython.display.Audio(f)